#!/bin/sh
#PBS -V
#PBS -j oe
#PBS -l nodes=1:ppn=1,mem=10gb,walltime=120:00:00
#PBS -q batch
Normal=$1
config_path=$2
source ${config_path}/config.sh

######################################################################################
## 原始的VCF文件
## 只考虑IM/IGC/DGC
all_tumor=`cat ${config_path}/Tumor_Normal_DNA.tsv | grep -w $Normal | \
awk -F'\t' '{print vcf_path"/"$3"_"$2"_GGA_Filter.vcf.gz"}' vcf_path=${vcf_path} | tr "\n" "\t"`

## 对每个vcf的多等位进行标化
## 拆分多等位
for file in `echo $all_tumor | tr '' '\n'`
do
out_file=`echo $file | awk -F'/' '{print $NF}' `
${bcftools} norm \
--fasta-ref ${seq_ref}/GRCh37.fa --multiallelics -both \
--output ${tmp_path}/${out_file}.forbcftools.vcf \
${file} 
## 压缩
${bgzip} -f ${tmp_path}/${out_file}.forbcftools.vcf
## 建立索引
${tabix} -f ${tmp_path}/${out_file}.forbcftools.vcf.gz
done

######################################################################################
## 标化以后的vcf
all_tumor_2=`cat ${config_path}/Tumor_Normal_DNA.tsv | grep -w $Normal | \
awk -F'\t' '{print tmp_path"/"$3"_"$2"_GGA_Filter.vcf.gz.forbcftools.vcf.gz"}' tmp_path=${tmp_path} | tr "\n" "\t"`

## 多等位标化
${bcftools} merge $all_tumor_2 > ${tmp_path}/Combine_${Normal}_GGA.vcf
#vcf-merge $all_tumor > ${tmp_path}/Combine_${Normal}_GGA.vcf
${bcftools} norm \
--fasta-ref ${seq_ref}/GRCh37.fa --multiallelics -both \
--output ${tmp_path}/${Normal}_GGA.forbcftools.vcf \
${tmp_path}/Combine_${Normal}_GGA.vcf 


######################################################################################
## 提取VCF对应的列，并且对没有call出来的突变赋值为0/1;100,0，保证所有突变都可以使用
cat ${tmp_path}/${Normal}_GGA.forbcftools.vcf | \
## add PASS
awk -F'\t' '{OFS="\t"}{if($0!~"#"){$7="PASS"};print}' | \
## 只取GT:AD:DP
awk -F'\t' '{OFS="\t"}
{
if($0~"#"){
print
}else{
split($9,spl,":");len=length(spl);
for(i=0;i<=NF;i++){
if(spl[i]=="GT"){GT=i};
if(spl[i]=="AD"){AD=i};
if(spl[i]=="DP"){DP=i};
}
$9="GT:AD:DP";
for(i=10;i<=NF;i++){
split($i,spl,":");
$i=spl[GT]":"spl[AD]":"spl[DP];
};
print
}

} 
' | \
## 非所有样本都能call出来的突变，补为0/1:100,0:100
awk -F'\t' '{OFS="\t"}{
if($0~"#"){print}
else{
str="";
for(i=10;i<=NF;i++)
{
if($i=="./.:.:." || $i=="././.:.:."){
$i="0/1:100,0:100";
}
};
if(str!="xx"){
print
}
}

}' | \
## 去除多等位非所有等位都能call出来的突变
grep -v "\.:" \
> ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf

######################################################################################
## 替换样本名
cd ${Tree_path}/
for sample in `cat ${config_path}/Tumor_Normal_DNA.tsv  | grep -w $Normal | awk -F'\t' '{print $5","$3}' `
do
echo $sample
Tumor=`echo $sample | awk -F, '{print $2}'`
Class=`echo $sample | awk -F, '{print $1}' | tr '-' '_'`
Class_new=$Class

line=`cat ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf | grep -n "CHROM" | awk -F":" '{print $1}'`

sed -i "${line}s/${Tumor}/${Class_new}/" ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf
done

Tree_mut_num=`cat ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf | grep -v "#" | wc -l`
Raw_mut_num=`zcat ${vcf_path}/Combine_${Normal}_forGGA.vcf.gz | grep -v "#" | wc -l`
echo ${Normal},$Raw_mut_num,$Tree_mut_num >> ${qc_path}/Treeomics_Mut_Num.list

## 检查突变是否一致
# ${bgzip} -c ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf > ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf.gz
# ${tabix} ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf.gz
# ${bcftools} isec -n+2 -c all -p ${tmp_path}/${Normal}_Tree ${Tree_path}/Combine_${Normal}_GGA_PASS.vcf.gz ${vcf_path}/Combine_${Normal}_forGGA.vcf.gz
